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Abstract 

Eigen's quasispecies system with explicit space and global regulation is considered. Limit 
behavior and stability of the system in a functional space under perturbations of a diffusion 
matrix with nonnegative spectrum are investigated. It is proven that if the diffusion matrix 
has only positive eigenvalues then the solutions of the distributed system converge to the 
equilibrium solution of the corresponding local dynamical system. These results imply that 
the error threshold does not change if the spatial interactions under the principle of global 
regulation are taken into account. 
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1 Introduction 

Molecular evolution is a complex nonlinear process that requires mathematical modeling to 
properly estimate mutual influence of the evolutionary forces such as mutation and selection. 
One of the simplest mathematical models describing the process of molecular evolution is due 
to Manfred Eigen [T7], see also [TJ [TBI 122] . It is worth noting from the very beginning that 
the classical models of molecular evolution, including the Eigen model, are usually written as 
systems of ordinary differential equations (ODE), and therefore describe mean field dynamics. 
It is well known that explicit spatial structure can significantly change the model behavior, see, 
e.g., [16]; therefore, it is important to evaluate the consequences of spatial variables. 

There are a number of studies that offer different approaches to molecular evolution models 
with spatially dependent interactions (see, e.g., [H [5j [121 EH CE3 OS E2]). Here we propose 
a general analytical framework along the lines we considered elsewhere (9[ [25], to tackle the 
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Eigen model with explicit spatial structure under global regulation principle. We start with the 
classical definition of the Eigen model. 

Consider a population of individuals (macromolecules) each of which belongs to one of k 
different types. Assume that the fitnesses (the replication rates) for these types are given by 
the vector a = (ai, . . . , and denote A the diagonal k x k matrix with a.^ % = 1, . . . , k on the 
main diagonal. Assume that upon reproduction for a macromolecule of type i it is possible to 
produce individuals of type j, and the probability of this event is qjf, hence qu = 1 — Ylj^ilji 
is the probability that no mutations occur at the reproduction. Let Q be the stochastic matrix 
with entries g^-. 

If n>i(t) denotes the absolute size of the population of the i-th type macromolecules, then, 
taking into account reproduction events with different fitnesses and mutations, assuming that 
the time is continuous in the model, after simple bookkeeping we obtain 

k 

hi = J ^2q ij a j n j , i = l,...,k, (1) 

3=1 

where the dot means the derivative with respect to the variable t. System (JTJ) can be written in 
a concise matrix form as 

h = QAn = Q a n, (2) 

where we put Q Q = QA and n(t) = (ni(t), . . . ,rik(t)). Note that for simplicity we suppress the 
dependence on t in the differential equations we consider. 

It is customary to consider not the absolute sizes of the corresponding types, but the fre- 
quencies 

nAt) 

Wi(t) = — k — , i = l,...,k, 

such that 

fc 

J>i(t) = l, (3) 

i=l 

and vector w(t) = (u>i(t), . . . , u>k(t)) belongs to the simples = {x 6 R fc : x > 0, Xa=i x i = 
For the new variables we obtain Eigen's quasispecies system [18] 



k 

w i = '^2 a jQij w j - Wif\t), i = l,...,k, (4) 

3=1 



with the initial conditions 

w(0) = £ = (&,...,&) €S fc . 

Here f (t) is the average fitness: 



A' 



f(t) = (a, w{t)) = mwi. (5) 
i=l 

Throughout the text we use (■, •) to denote the inner product in R fc . 
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It is possible to rewrite system in the form of an almost linear system if we add the 
following nonlinear equation for the average fitness 



where (H|) was used. 

The only rest point of dynamical system ^) can be found as the solution of the following 
problem 



where I is the identity k x k matrix. 

Due to the fact that the matrix Q a can be assumed to be primitive, then, invoking the 
Perron-Frobenius theorem, we have that the unique solution of ([7]) is given by the positive 
eigenvector w = (wi, . . . ,wp.) € of Q a that corresponds to the dominant positive eigenvalue 
y = (w,a). Note that f l is also the rest point of ([6]), because ([6]) yields that (Q a w,a) = 



The two major findings of Eigen, who analyzed model (j4]), were that the evolution works 
not on separate types of macromolecules, but on the vectors composed of the mutant clouds; in 
evolutionary process the winner is not the type with the maximal fitness, as it would happen if 
all qij = 0, but the eigenvector w corresponding to the dominant eigenvalue. This vector was 
called the quasispecies, hence the name of the model. The second important observation is the 
presence of so-called error threshold (for more details see, e.g., [32]), which predicts that for many 
fitness matrices A for the replication rates there exists a critical mutation rate over which the 
distribution of frequencies of different types of macromolecules becomes uniform. That implies 
that no statistical method can determine which type has the maximal fitness (for a discussion see, 
e.g., [H O El [TOl E21 EJl EHl [32] ) . The discussion on the nature, exact definition and importance 
of the error threshold for the real populations is ongoing, and we do not plan to touch on these 
issues here. Instead, we concentrate on the following problem: How to incorporate a spatial 
structure into the model ([U), and what consequences can be expected if the spatial interactions 
are modeled. 

The rest of the paper is organized as follows. In the next section we discuss our approach to 
add space to model ([!]) and state, in a nontechnical language, the main results of our analysis. 
In Section [3] we specify the mathematical problem, the functional spaces in which we are looking 
for a solution, and discuss in which sense we use the word solution. Sections [4] and [5] provides 
the details of the proof of the main theorem, including some more subtle corollaries. In Section 
[6] we illustrate our findings numerically. 

2 Eigen's quasispecies model with diffusion. Main theorem 

It is a usual approach in mathematical modeling (e.g., [HIES]) to add the spatial structure to 
a system of ODE 




(6) 



(Qa ~ If(t))w = 0, 



(7) 



f\w,a) = (f l f. 




(8) 





rn 



(9) 
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where dtn = Qjjr, D is the matrix of diffusion coefficients, m = 1,2, or 3 depending on the 
dimension of the considered physical space (line, plane or space) and A is the Laplace operator. 
System ([9]) should be supplemented with appropriate initial and boundary conditions. It is well 
known that the diffusion can be a destabilizing force in a sense that an asymptotically stable 
rest point of ([8]) can become an unstable spatially uniform rest point of ([9]) through the Turing 
bifurcation |19t 129] . The exact conditions whether the Turing bifurcation occurs depend on a 
particular system. Therefore, we are interested in a possibility of a Turing bifurcation for a 
spatial analogue to system (jlj). 

An important observation here is that the general approach that leads from ([8]) to ([9]), 
however, will not work for the quasispecies model (j3J), since we have a second additional condition 
that w(t) € Sf. for any time moment t; the simplex Sk is invariant in R fc . 

There are several possible approaches to circumvent the mentioned obstacle and to add 
explicit spatial structure to (jU , for a short review see [25J . We believe that the correct approach 
is to follow the lines how the system for the absolute sizes ([TJ was transformed into the system 
for the frequencies (jl]) , (this correspond to the general ideology of the theory of selection systems 
pi]). Using this approach we arrive (for the details see Section [3]) at the equation 

k 

dtUi = ^2 a jQij u j ~ u if s (t) + dijAuj, x G $7 C R m , i = 1, . . . , k (10) 
j=i j 

for the new variable u(x, t) = (u%(x, t), . . . , Uk(x, t)). In (fTUj) we use notation 

/'(*) = E / ^ u ) dx ( n ) 

for the mean integral fitness that guarantees that in the domain the principle of global regu- 
lation 



V] / Ui(x,t) 



dx = 1 



holds for any time moment t. 

This approach was used in [TJ [9] for various classes of replicator equations. For the Eigen 
model ([4]) this approach in the case when space has one dimension was applied for the first 
time in [31], and it was argued that the spatially uniform steady state for the local model (j3|) 
remains asymptotically stable for the case of system (jlOp when the diffusion matrix is diagonal. 
We generalize this statement for an arbitrary dimension of the physical space and any diffusion 
matrix D that has positive eigenvalues. The main statement is given below, and the rest of the 
paper is devoted to the proof of it. 

Theorem 1. Consider simultaneously the local Eigen quasispecies model and the spatially 
explicit model (llOp with zero Neumann boundary conditions (the boundary ofQ is impenetrable). 
Let domain Vt be bounded and have a piecewise smooth boundary, diffusion matrix D have only 
distinct positive eigenvalues. Then 

1. There exits only one steady state solution to (jlpp . and it coincides with the rest point w 
of ®; 
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2. This steady state solution w to (jlOp is globally stable. 

Remark 2. The convergence for the solutions of system (|1U|) to the rest point w is understood 
as the convergence in an appropriate functional space, the details are presented in Sections [3] 
andU 

Remark 3. Model ([!]) is not the only possible mathematical reincarnation of the Eigen qua- 
sispecies model. Another, and somewhat more justified approach to deduce a model that is 
described by a system of ODE is to assume that the mutations and reproduction events are 
separated, and reproduction proceeds error-free. In this case, it is possible to show pQ that the 
system takes the form 

k 

m = {mi - f(t))wj + fajWj, i = l,...,k, (12) 

i=l 

where now mi are Malthusian fitnesses (the reproduction rates) and faj are the mutation rates, 
Mii = ~~ Mj*> an< ^ the mean fitness is given by 

k 

/'(*) = ^rniWi. 

i=l 

Systems (jj]) and (|12p have very close properties [21j . Theorem [T] (with necessary modifications) 
remains valid for system (fT2|) as well. 

3 Notation and derivation of the Eigen model with diffusion 

In this section we present a general approach to derive a spatially explicit counterpart of the 
classical Eigen's quasispecies model 

Let C R m , m = 1, 2, 3 be a bounded domain with a piecewise-smooth boundary 7. Denote 
n(x, t) vector-function that describes absolute sizes of different macromolecules composing the 
system; in coordinates n(x,t) = (ni(x,t), . . . ,rik(x,t)), t > 0, x G $7 in space R m . 

Consider matrix D = (dij)kxk with positive eigenvalues fa > 0, i = 1, . . . , k and assume that 
there exists a matrix U such that 

D = UMU-\ (13) 

where M is a diagonal matrix with fa on the main diagonal. The diffusion matrix D describes 
the action of spatial diffusion on macromolecules; note that we do not require matrix D be 
diagonal, and hence include in consideration possible cross diffusion interactions that describe 
a directed movement of macromolecules with respect to an external signal (which can be, e.g., 
chemical, see for more detail [51 [20]). 

Following the standard approach (e.g., |TT|) the dynamics of the populations of macro- 
molecules can be described with help of the following PDE 

d t rii(x,t) = (Q a n(x,t))i + (DAn(x,t))i, i = l,...,k. (14) 
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Here A is the Laplace operator, in Cartesian coordinates 

i=l 1 

Q a is the same matrix as in equality (|2|), and the following notations are used in the rest of the 
paper: 

{Q a n(x,t))i = ^2^qijajnj(x,t), 
u 

(Mn(i,t))i = ^ ^ dijArij(x, t). 

The initial data for the concentrations of macromolecules are given by 

ni(x,0) = i = 1, ...,k. 

It is natural to assume that we consider a closed system, i.e., we have zero flow boundary 
condition (the Neumann boundary conditions) 



drii(x, t) 



= 0, i = l,...,k, 

xS7 



where v is an outward normal to the boundary 7. 

We would like to rewrite system (|14p for the relative concentrations. Consider 



Ui(x,t) = — r , 1 = 1, ...,k. (15) 

Ej=iln n j( x ^) dx 

We obtain equality ^ in the form of the following integral invariant 

n „ 

i(u, i k )) = ^2 Ui<yX > ^ dx = lj ( 16 ) 

where the following notation 

((«, l fc )> = f (n(x, t), l fc ) dx, l fe = (1, . . . , 1) G R fe 
Jn 

was used; in words we use the angle brackets {F) to denote the integral of F over Q, and the 
parentheses, or round brackets (•, •) to denote the inner product in R fc . 
From (1151) it follows that 



d t ni((n, l fe )) - ni((d t n, l fc )) 
dtUi = 77 — T^2 = 

_ ((Q a n)i + (-DAn)j) ((n, l fc )> - n;((Q a n + OAn, l fe )> 
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After simplification, we obtain the following equations 



dm = (Q a u)i - Uif s (t) + (DAu) h i = l,...,k, (17) 
with the initial and boundary conditions 

dui{x, t) 



Ui(x,0) = 4>i(x), 



= 0, t>0,i = l,...,k. (18) 



dv 

Function f s (i) is the average integral fitness of the system and is given explicitly by 

f s (t) = {(u(x,t),a)). (19) 
For obtaining the last expression we used boundary conditions (|18p and Green's formula 



f 9ui(x, t) 

Aui(x,t)dx = / as = 0, % = 1, . . . , k. 

J 1 dv 

Suppose that for any fixed t each function Ui(x,t) is differentiable with respect to t and 
belongs to the Sobolev space W 2,m (£l) if m = 1, 2 or W 2,2 (Q) if m = 3 as the function of x for 
any fixed t > . Here W m > 2 , m = 1,2 is the space of functions, which have square integrable 
derivatives with respect to variables Xi up to the order m. We require this condition because, 
due to the embedding theorem we have that any function from the space W s ' 2 (£l) is continuous, 
except possibly on the set of measure zero (see, e.g., [191 for details). 

Denote £lt = ^ x [0, +oo) and consider the space B(Qt) of functions with the norm 

\\v\\ B = max {11^11^,2 + 119^11^/^,2} 

Denote Sk (^t) the set of non-negative vector-functions u(x,t), such that Ui(x,t) £ B (Q t ) , i = 
l,...,k that satisfy (fT6|) . In the following we look for weak solution to (fT7|) . (fTKj) . i.e., the 
solutions that satisfy the integral identity 

poo r f'OO r 

dtiiir)(x,t)dxdt = / / (Q a u(x,t))ii](x,t)dxdt — / / (DVu,Vt])idxdt 
Jo Jn Jo Jn 

for any function r](x,t) on compact support in variable t, which is differentiable on [0, +oo) with 
respect to t and belongs to the Sobolev space W S ' 2 (Q.) for any fixed t > 0, s = 1, 2. 

We remark that system (|17p is not a usual system of PDE because its right hand side contains 
the functional f s (t) on solutions Ui(x,t). Nevertheless, it is possible to consider this system in 
the classical sense if we add the following ODE for the average integral fitness that can be 
obtained similarly to equation ([6]) 

^(t) = ((Q a u,a))-(nt)) 2 . (20) 

The steady state solution to (|17p can be found as a solution to the following PDE of elliptic 
type: 

{DAv{x)) i + {Q a v{x)) i -vi(x)f = 0, i = !,...,*, (21) 
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with the boundary conditions 

dviix) 



dv 

The integral invariant (|16p now reads 



= 0, i = l,...,k. (22) 

££7 



((v(x), l k )) = V f Vi (x)dx = 1. (23) 
i=i 7n 

Denote by Sfc(fi) the set of all non-negative vector function v(x) = (vi(x), . . . ,Vk(x)), Vi{x) G 
W s ' 2 (n) (s = 1 or s = 2) that satisfy fl23]). Using {TS} and ([23]) we obtain that 

/- = ((«(x),a)), 

i.e., f s is a constant. 

4 Steady state solution to the distributed Eigen's system 

In this section we prove the first part of Theorem [TJ 

Theorem 4. If the diffusion matrix D has only positive eigenvalues \n, i = l,...,n, then the 
steady state solution to the system (j21j) . (I22p in Sk (fi) coincides with the rest point w of the 
local system 

Proof. Consider the following boundary eigenvalue problem 

dip 



Atp(x) = -Xip(x), 



= 0. (24) 

XG7 



The eigenfunction system is given by ipo{x) = 1, {ipj(x)}JL 1 and forms a complete system in the 
Sobolev space W s ' 2 (Q); moreover, we have 

ipi(x)ipj(x) dx = Sij, (25) 

n 

where 5ij is the Kronecker symbol. The corresponding eigenvalues are known to satisfy 
Ao = < Ai < A2 < ... < \j < ■ ■ ■ , lim Xj = +00. 

Consider solution v(x) = (v\(x), . . . ,Vk{x)) G 5^(0) to the system (|2ip. ([22^1 . Assume that this 
solution has the following form 

00 

Vi{x) = 4 + ^4^{x), i = l,...,k (26) 
i=i 

which is always possible since the eigenfunctions {t/^j (x)}?2 =0 of (|24p form a complete system in 
W S,2 (Q). Here c?, cj are constants, i = 1, . . . , k, j = 1, 2, 
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We substitute (j26|) into (f2Tj) and multiply consecutively the resulting expression by the 
eigenfunctions of (|24p . After integrating over Q we obtain a systems of linear equations for 

(QaC°)i - fcf = 0, i = l,...,fc, (27) 

(QaC^i - /'cj - (DA^-C^i = 0, i = 1,2,... (28) 

where Xj are the eigenvalues of ([Ml), C° = (c§, . . . , eg) and C J ' = (cj, . . . , c£) are vectors in R , 
j = 1,2, . . . Here we essentially used the orthogonality of the eigenfunctions ipj(x). 
The integral invariant (|23p and the average integral fitness take the form 

i=l i=l 

Equation (|27|) has the same form as the equation for rest point of system 0. Therefore, the 
unique solution of that equation coincides with the coordinates of rest point, i.e. C° = w, which 
implies 

r = f. 

Now consider system ()28|) . Since D satisfies (|13p . then there exists a non-singular matrix U 
such that U~ 1 DU = M, where M is a diagonal matrix with the entries /ii, . . . , /i& on the main 
diagonal. Consider the transformation 

& = UY\ j = 1,2, ... 

where = (y-[, ■■■,yV) is a vector in R fc . Equations (|28|) take the following form 

RoX^ = (f l I + \jM)Y-i , j = 1,2,... (29) 

The eigenvalues of the matrix 72 a = U~ 1 Q a U and matrix Q« coincide. 

Using the fact that /' is the maximal eigenvalue of the matrix Q a we obtain 

((C/^QaC/)^',^') < f l \\Yi\\ 2 . 

By the assumptions, the eigenvalues fj,{ of the diffusion matrix D are positive and hence 

{(f l I + \ j M)Yi,Yi)>f\\Yi\\ 2 . 

Therefore, equation (|29|) has no nontrivial solutions. It implies that the only solution to equation 
(|27l) is exactly the rest point of the local system ([4]) . □ 

Remark 5. The case when at least one of the eigenvalues of D is zero is more intricate. Assume, 
without loss of generality, that = 0, then it is possible that 

(MYi,Yi) = Q, i = 1,2, ... . 

Therefore, the equation ([29]) can have a countable set of solutions. 
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Consider the following example. Suppose that the rest point of @ is an eigenvector of D 
with the zero eigenvalue: Dw = 0. Then the vector functions 

1 

101 



^) = TF^( 1 + ^(z)), i = 1,2, . . . 



are the solutions to the system (|2ip . (|22p . Here \Q\ is the measure of f2. Indeed, we have 

DAv j + Q a v j - If l v j = —(Dw)ijji(x) = 0, 



— / ^(x)da; = — / T&ids 



5 Stability of the spatially homogeneous solution to the dis- 
tributed Eigen's system 

In this section we provide a proof to the second part of Theorem [TJ 

For the following we will need the following fact: If the system x = Ax{t) with a constant 
matrix A has bounded solutions when t — > oo (i.e., the origin is stable) then the system y = 
(A + B(t))y(t) also has bounded solutions if J* °° ||-B(i)||<ii < oo; the proof can be found in, e.g., 

Theorem 6. // the eigenvalues of the diffusion matrix D are positive then for any initial data 
4>i(x) £ S n ($7) solutions to the system ()17p . f)18|) converge in S n (£lt) to the solutions of the 
system Q with the initial data 

£i= 4>i(x)dx, i = l,...,k. (30) 
Jn 

Proof. As in the proof of Theorem [4] we look for the solution Ui(x,t) G S n (Slj) of (fT7l) . (fT8j) in 
the form of the following expansion 

oo 

u i (x,t) = c i (t) + ^24(t)^ J (x), i = l,...,k. (31) 

3=1 

Here cj (t) are functions of t, i = 1, . . . ,k, j = 1,2,.... This expansion is possible because the 
system of eigenfunctions {ipj(x)}j^ forms a complete system in W S,2 {Q), s = 1,2. We obtain 
the following system of ODE for the unknown (t) 

^(t) = (Q*c (t)) l -f s (t)c°(t), (32) 
dd ■ 

= - nt)4(t) - (Dx^m, (33) 

where A,- are the eigenvalues of ((Ml); c°(t) = (c?(t), . . . , c° k (t)), c j (t) = (c[(t), . . . , d(t)). The 
integral fitness (|19p has the form 



i=i 



10 



Since 



4>i(x)dx = / Ui(x,0)dx = & 
Jn 



then = fa, i = 1, ... ,k. From the uniqueness theorem it follows that the solution to the 
systems (jl]) and (i32j) coincide, i.e. c®(t) = W{(t) and therefore 

/*(*) = fit)- (34) 

Now consider system (|33|) . As in the proof of Theorem [J] consider the transformation U which 
reduces the diffusion matrix to the diagonal form. System (|33|) takes the form 

= {U~ l Q a U)Yi(t) - (f l (t)I + XjM)Y^(t), j = 1,2,... 
where M is the diagonal matrix. Using 

f\t) = {f\t)-j l ) + J\ 

we obtain 

^ = AY\t) + B{t)Y\t), 

where A is the constant matrix 

A=(U- 1 Q a U)-(f s I + X 1 M), 
B{t) is a diagonal matrix with the elements 

hi{t) = f - f\t), i = l,...,k. 
Recall that f(t) is given by ([5]) and f l = (w,a). To investigate the system's behavior for 

— =AX j (t), X j (t)eR k 
we introduce the following Lyapunov functions 

V(X*) = ±(X3,X3) = l\\Xi\\ 2 . 

One has 

V(X>) = (AX j , X j ) = ((U^QatyX^X') - ((f l I + \jM)Xi,Xi). 
On the other hand, 

((U- 1 Q a U)X j ,X s ) < f\\X j \\ 2 , 
((f l I + X j M)X j ,X j ) > /'||X J '|| 2 . 

Therefore the equilibrium X 3 = is asymptotically and exponentially stable. 
Now estimate the following integral 

oo 

\T-f{t)\dt. (35) 

o 
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It can be reduced to 

/•oo 

/ \wi - Wi(t)\dt, i = l,...,k, 
Jo 

where Wi(t) are the coordinates of the solution to the equation and u)j are the coordinates 
of the rest point of (jU). 

Using the well known results for the convergence of the solutions to the rest point of (jJJ) 
(e.g., |23j) we obtain that there exist positive constants K and k such that 



\wi — Wi(t)\ < Ke 



-Kt 



Therefore the integral (|35p converges. It implies that the solutions to the system (j33|) tend to 
zero when t — > oo. This completes the proof. □ 

Corollary 7. If at least one eigenvalue of diffusion matrix is equal to zero then it is possible to 
have a countable set of solution u 3 (x,t) £ S n ({lt), j = 1,2,... to the system ([17p . (|18p . These 
solutions converge to the solutions of the dynamical system @ in the sense of average integral 
value, i.e., 

lim / u J Ax,t)dx — ?> Wi(t), i = l,...,k, j = 1,2, ... 
Jn 

where w(t) is the solution of the system @. 

Indeed in this case we can not guarantee convergence the solutions of (|33[) to zero. Consider 
a countable set of vector functions 



u?(x,t) = -r^-(u(t) +U)lf)j(x)), 



DU j (x,t) + Q a U j {x,t) - If l u j (x,t) - d t u j = 
T^-((Q a uj(t) - lfoj(t) - d t w) + (D6j + Q a LO - If l u)i>j{x)) = 0, 



since DCj = 0. Additionally. 

k 



[ uj(x,t)dx = -^ ^T j ui(t)dx = ^Wi(t) = 1. 

i=l Jn ' ' i=l i=l 



6 Numerical calculations 

We illustrate the results of analytical calculations with the following example. Consider the 
population of macromolecules that have length L = 4 and composed of ones and zeros. As- 
suming that the corresponding fitnesses depend only on the amount of zeroes and not on their 
positions (the fitness landscape is permutation invariant), the mutations are independent of 
the site position and occur with probability q per individual per position per replication, then 
the whole quasispecies system can be rewritten as L + 1 dimensional system with the variables 
u(x,t) = (u x , . . . ,u 5 ). 

Numerical solutions for such system are shown in Fig. [TJ In particular, the top panel gives 
the time evolution of the spatially independent frequencies Ui(t) = f n Ui(x, t)dt, and the bottom 
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panels show how spatially non-homogeneous initial conditions becomes homogeneous under the 
influence of spatial diffusion. The overall picture is in agreement with the analytical investigation 
from the previous sections; the solutions approach the values w of the rest point of the local 
Eigen system, which can be found as the coordinate of the positive eigenvector. 
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